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Abstract 

It is known that the exact analytic solutions of wave scattering by a circular cylin- 
der, when they exist, are not in a closed form but in infinite series which converges 
slowly for high frequency waves. In this paper, we present a fast numerical solution 
for the scattering problem in which the Boundary Integral Equations, reformulated 
from the Helmholtz equation, are solved using a Fourier spectral method. It is shown 
that the special geometry considered here allows the implementation of the spectral 
method to be simple and very efficient. The present method differs from previous ap- 
proaches in that the singularities of the integral kernels are removed and dealt with 
accurately. The proposed method preserves the spectral accuracy and is shown to have 
an exponential rate of convergence. Aspects of efficient implementation using FFT are 
discussed. Moreover, the boundary integral equations of combined single and double- 
layer representation are used in the present paper. This ensures the uniqueness of the 
numerical solution for the scattering problem at all frequencies. Although a strongly 
singular kernel is encountered for the Neumann boundary conditions, we show that the 
hypersingularity can be handled easily in the spectral method. Numerical examples 
that demonstrate the validity of the method are also presented. 


This work was supported by the National Aeronautics and Space Administration under NASA Contract 
NAS1-19480 while the author was in residence at the Institute for Computer Applications in Science and 
Engineering, NASA Langley Research Center, Hampton, VA 23665, USA. 
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INTRODUCTION 

The exact analytic solutions of wave scattering by a circular cylinder, obtainable 
for simple incident waves, are not in a closed form but in infinite series of Bessel and 
Hankel functions of increasing orders. Such solutions converge slowly, especially for 
high frequency waves, which render their numerical evaluation inefficient. This paper 
presents a fast numerical solution of wave scattering that only requires the computation 
of Bessel and Hankel functions of order zero. Furthermore, the numerical solution is 
valid for any form of the incident waves of all frequencies. 

When developing numerical solutions, wave scattering problems are often conve- 
niently formulated in Boundary Integral Equations (BIE)*. The advantages of the 
Boundary Integral Equation Method (BIEM) include reducing the dimension of the 
problem and transforming an infinite domain to finite boundaries in which the far field 
radiation condition is satisfied automatically. The Boundary Integral Equations are 
commonly solved computationally by the Boundary Element Methods (BEM) 2 . In 
this method the boundary is divided into finite elements and integrations over each 
boundary element are approximated by quadratures, e.g. the linear elements. 

In this paper, we develop a Spectral Method of solving the Boundary Integral 
Equations, reformulated from the Helmholtz equation, for numerical solutions of wave 
scattering by a circular cylinder. Previously, for this special geometry, a “fast numerical 
method” based on the Fourier approximations has been formulated by Bojarski 3 , who 
pointed out that the boundary integral equation of wave scattering can be solved easily 
and efficiently in the Fourier spectrum domain of the solution. Due to the simplicity 
of the geometry, an explicit relation between the Fourier coefficients of the solution 
and those of the boundary condition was found. It was argued that the numerical 
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approach was more efficient than directly evaluating the infinite series of the exact 
solutions. Indeed, the exact solutions contain Bessel and Hankel functions of higher 
orders whose numerical evaluation is more difficult and costly as the order increases. 
Recently a similar approach has been used and extended by Schuster 4 for a wave 
transmission problem of concentric cylinders. 

In the present paper, we point out that the numerical formulations given previously 
are not achieving the optimal accuracy of the Fourier spectral methods. It is known 
that, although any periodic function can be approximated by a truncated Fourier se- 
ries, the rate of convergence of such an approximation depends on its smoothness. 
Unfortunately, the integral kernels for the Helmholtz equation are not smooth. In par- 
ticular, the 2-D Green’s function of the Helmholtz equation, appearing in the integral 
equations, possesses a logarithmic singularity. Furthermore, the normal derivative of 
the Green’s function also contains a term involving the logarithmic function. The non- 
smoothness of the integral kernels, however, was not explicitly treated in the previous 
formulations. It will be seen that it is critical to remove the non-smoothness of the 
integral kernels in order to achieve fast convergence in the Fourier spectral formulation. 
By a proper treatment of the singularities, the present numerical formulation yields 
accurate solutions with significantly fewer datum points. Moreover, the boundary in- 
tegral equations of combined single and double-layer representation are used in the 
present paper. This ensures the uniqueness of the numerical solution for the scatter- 
ing problem at all frequencies^. Although a combined layer formulation results in 
a strongly singular kernel for the Neumann boundary conditions, we show that the 
hypersingularity is handled easily in the spectral method. 

In the next section, the formulations of the Boundary Integral Equations for wave 
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scattering problems are given. Then, in sections II and III, the Fourier spectral methods 
for the Dirichlet and Neumann boundary conditions are presented. Numerical results 
are shown in section IV. Section V contains the conclusions. Some analytic results are 
also given in the Appendix. 


I. BOUNDARY INTEGRAL EQUATIONS 

Let us consider wave scattering by a circular cylinder T of radius a. The wave 
equation for the scattered function <}> with assumed time dependency of e~ twt is reduced 
to the Helmholtz equation 

V 2 <t> + K 2 <j) = 0 (1) 


where k 

d 2 d 2 

dx 2 + dy 2 


= u/c ( c is the wave speed) and V 2 is the 2-D Laplace operator V 2 = 
. The boundary condition considered in this paper will be one of the following 


types : 


Dirichlet(soft) : <f)(r ) = 6(r ) On T 


or 


Neumann(hard) : — (r ) — b(r ) On T 

The Helmholtz equation (1) together with the boundary condition can be refor- 
mulated into a Boundary Integral Equation. This can be done in various ways 1 ’ 5 . For 
scattering problems considered in the present paper, we use a combination of single 
and double-layer formulation in which the solution <f> at any point r 1 in the scattered 
field is represented by an integral on the boundary as 5 


H?') = jf /(r-)<JT 


( 2 ) 
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where 77 is any real number such that 


T)Re{n ) > 0 

The use of a combined formulation ensures the uniqueness of the numerical solution 
for exterior problems 1,5 . In (2), /(r) is an unknown layer distribution function and 
the Green’s function G(r, r'), whose form will be given later, satisfies the following 
equation 

V 2 G + K 2 G=-6{r-r') (3) 

Q 

Here the normal derivative -7— is assumed to be taken in the direction outward from 

an 

the cylinder. 

The Boundary Integral Equation associated with the layer representation (2) is 5 

\l(0 + l(^~ w) m cW = HO (4a) 

for Dirichlet boundary conditions and 

7'«> + 1 GH - *5?) WMT = Mr;) (46) 

for Neumann boundary conditions , respectively. In (4a) and (4b), r r denotes the bound- 
ary points. After the layer distribution function / has been solved from the integral 
equation (4a) or (4b), the solution of the Helmholtz equation ^ is found by the bound- 
ary integral (2). 

Now for a circular cylinder of radius a, the boundary contour can be expressed as 
f r (0) = (acos0, asinfl), 0 < 0 < 27t (5) 

The normal vector to be used in (4a) and (4b) is n = (cos 9, sin 0). 
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The Green’s function and its normal derivative are 5,6 


and 







dG 

dn 




( fr( 0 )-r r (6 ')) n 
|r r W-f r (^)| 


( 6 ) 


= [ 2«a 

4 1 


e-e’ 


sin 


e-e' 


sin 


(7) 


sin 


0 - 0 * 


in which we have used the fact that |r r (0) - r r (0 t ) \ = 2 a 

It is important to note here that G and are functions of 9 — 0 f . As will be seen 
later, this allows the implementation of the Fourier spectral method to take a simple 
form. 

We thus express the boundary integral equation (4a) for the Dirichlet boundary 
conditions as 


1 f 2rr f dG \ 

2 /(*') + J - tf) - ir)G(0 - 0')) f(9)ad9 = b(O') (8a) 

and equation (4b) for the Neumann boundary conditions as 

+ 1 ( 6 - ^ - ‘1^ - <’')) /(*)« de = W) («&) 

For clarity, the dependencies on 9 and 9 1 have been expressed explicitly in (8a) and 
(8b). 

In the next two sections, we give the numerical formulations of solving the inte- 
gral equations (8a) and (8b) by a Fourier spectral method. Since different types of 
singularities are encountered, the two equations will be dealt with separately. 
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II. SPECTRAL METHOD FOR DIRICHLET BOUNDARY CONDITIONS 


A. Formulation 

Let the layer distribution function f(9) and the boundary condition b(0) be ap- 
proximated by the truncated Fourier series as 

N/ 2-1 

m = E (9) 

n=-N/2 
N/ 2-1 

b(e)= Ke ' n9 ( 1 °) 

n=—N/2 

where b n are obtained by the FFT from prescribed boundary condition and /„ are the 
unknown coefficients. In (9) and (10), the particular form of truncated Fourier series 


has been taken for the convenience of applying FFT programs. 

Substituting (9) and (10) into the boundary integral equation for the Dirichlet 
boundary conditions (8a), we get 

Nf 2-1 N/ 2-1 f 2* /f)r \ I 

l E /.«** + E U/ (^;(e-e')-iva(e-e'))e^ade = E 

" mr/o L ^0 \ / J n=— jV/2 


n=-N/2 


n=-N / 2 


(ii) 


For simplicity, let 


x = 9 — 9 ! 


By equating the coefficients of e tn ^ , equation (11) is easily reduced to 

1 fn + fnj^ (^(x)-iriG(x)^e ,nx adx = b n (12) 

for —N/2 <n<N/ 2-1. 

It is seen that the integral appearing in (12) are related to the Fourier coefficients 
°f |g(x) and ^(a;). From (6) and (7), it is also clear that both are periodic functions of 
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x , with a period of 2i r. Thus if we let G{x ) and ^(x) be approximated by truncated 


Fourier series as 


N/ 2-1 

G(x)= Y. 

n=—N/2 


dC 

s^*)- E 

n=—N/2 


then, the integral in (12) equals to 2Tra(h n — ir)g n ). It follows that 


(13) 


(14) 


-/„ + 2naf n (h n - irjg n ) = &„ 


(15) 


Therefore, the Fourier coefficients of the layer distribution function /(<?) are ob- 
tained explicitly as 


In — ] 


bn 


k + 2na(h n - irjg n ) 


(16) 


The above equation shows that once the Fourier coefficients of (T(x) and f^(^) 
have been found, the layer distribution function f(0) is known immediately. 

Actually, the Fourier coefficients of G(x) and §^(z) can be found in exact form 
using higher order Bessel and Hankel functions. They are derived in Appendix A. 
Nonetheless, the numerical evaluation of the exact expressions becomes more ineffective 
and costly as the order of the special functions increases. In what follows we give the 
numerical method that computes the Fourier coefficients g n and h n accurately and 
efficiently. 


B. Computation of g n and h n 

In general, the Fourier coefficients of a periodic function can be obtained efficiently 
by using a Fast Fourier Transform algorithm (FFT). However, the accuracy of the 


7 



Fourier coefficients computed by the FFT using a given number of datum points de- 
pends on the smoothness of the function. Only when the function is infinitely smooth 
(i.e. infinitely differentiable), the error of Fourier coefficients computed by FFT de- 
cays faster than any power of 1/N, where N is the number of datum points. Such a 
convergence is often referred to as an exponential convergence and the method is said 
to have spectral accuracy 7,8 . Our aim here is to compute g n and h n by the FFT with 
spectral accuracy even though the functions G and ^ are not smooth. 

In the numerical approaches proposed previously 3,4 , the Fourier coefficients g n and 
h n were computed directly as the FFT of the G(x ) and |^(x) respectively. However, 
the Green’s function G(x) has a logarithmic singularity at x = 0, where 6 = 6' , due 
to the Hankel function of order zero in (6), and its Fourier series converges at the rate 
of 1 /N. Thus direct computation of g n from G(x) using FFT yields results whose 
accuracy is only comparable to a first order method. Furthermore, the function |^(x) 
also has a non-smooth derivative at x = 0, and its Fourier series converges at the 
rate of 1 /N 3 . Thus direct computation of h n from x ) is only comparable to a 
third order method. Alternatively, as will be shown below, by properly treating the 
non-smoothness of G(x) and j£{x), g n and h n are computed with spectral accuracy. 

To examine the singularity of G(x), we note that 

G(x) — -H { 0 1) ^2/ca sin— h = — |jo ^2/ca sin— + iYq ^2/c a sin ^ ^ 


in which Jo and Yo are the zeroth order Bessel functions of the first and second kind, 

respectively. Using the asymptotic series for small arguments, we have 9 

z 2 z 4 

^W- 1 - 7 + 64-' 
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It follows that, for |x| small, 


G(x) = --Gn (/c a [sin | ) J 0 (2 ku |sin j|) _ ^ + J + 0(x 2 ) 


(17) 


in which 0(x 2 ) represents a power series in x 2 , and 7 is the Euler’s constant, 7 = 
0.577215... . To compute the Fourier coefficients of G(x) efficiently and accurately, we 
note that the Fourier series of the logarithmic periodic function In («a |sin ||) in (17) 


is 


6 . 


i" (“hf|)= ta (f)-E 


cos(nx) 


n=l 


n 


(18) 


Thus we can “subtract out” the singularity in G(x) by forming 


G(x) = — (2*0 sin — ^ — In (/ca |sin — Jo ^2/ca sin ^ |) (19) 


and then writing the Green’s function as 


G(x) = G(x) — — In (tea |sin — |) Jo (2 na sin ^|) 


(19') 


It is easy to see that G(x) is finite for all values of x. Furthermore, both G(x) and 
Jo (2/ca | sin |-|) in (19') are periodic and infinitely differentiable. Thus their Fourier 
coefficients can be computed with spectral accuracy using FFT. The Fourier coefficients 
of the Green’s function G(x), g n , will be computed according to (19') where the term 
involving the logarithmic function is computed by using convolution sums. 

We now study the non-smoothness of the normal derivative of the Green’s function 

dG 

-^•(x). The asymptotic series of the Bessel functions of first order for small argument 


are 


9 . 


T , . 2 2 

^ ~ 2 ” 16 + 


y , w = __ + _ ln (|) JlW + 


27- I 

2n 


-z — 


9 



Then 


dG . . in (i) ( . x|\ I . xl 

= -jH\ ( 2 *a s,n 2 1) j sm 2| 

= — T i*^ 1 (? Ka | s * n 7 1) * 1 ^ (^ Ka s ' n ^ 


. x 
sin — 
2 


+ 


£-\n ( 
2 ^ r V 


( . X 

V _ { . X \ 

. x 

[ kcl sin — 

J i 2 kcl sin — 

sin - 

V 2 , 

1 \ 2 J 

2 


+ 0(x 2 ) (20) 

4ira ' 2tt V ” I 21/ ~ * V I 2 1/ I 21 v ' v ' 
Thus although ^ is a finite function, due to the logarithmic function appearing in 
the second term shown in (20), it does not have a smooth second derivative at x = 0. 
For this reason, its Fourier approximation will converge only at the rate of 1/N 3 . 

The Fourier coefficients of however, can be found easily using the relation to 
g n given in Appendix A. In particular, we have 

'(s'n+l — 9n— l) n ^ 0, 2 ) "2 ^ 


/ K 2 a , 


K= { 


4n 

K 2 a 1 

—r{92 - go) - -g\ 

Q d 

n 2 a 


4ra g ~T +l 
n 2 a 


n = 0 


n = -f 


n = f-l 


( 21 ) 


— 9^-2 

An * 

Thus it is only necessary to compute g n , the Fourier coefficients of G(x). 


C. Fast Fourier Transforms 

The numerical implementation of computing g n by (19 ; ) is given in this subsection. 
Let us introduce Fourier collocation points 

j = 0,1,2 JV- 1 

For convenience of discussion, denote the following Fourier series approximations 

N/ 2-1 

G(x)= Y, 9ne~ inx (22 a) 

n=-N/2 
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N/ 2-1 

f I X l\ * 

Jo [Zita | sin - J = 2J Pne - *" 

n=-N/2 


( 226 ) 


The coefficients of these expansions are computed by FFX (backward in the usual 
sense) as follows: 


1 N ~'- 
s " = n'E 

j = 0 

Pn = ]jYl J ° ( 2 «a|s in y|) 

j = o 


e mx * 


(22a') 

(226') 


in which G(xj) is computed by (19). For the value of G(x) at x = 0, the following 
limit, obtained from (17), can be used : 

^( 0 )«-JL + i 

W 2tt 4 

In addition, we denote (18) as 

0 ° 

In \Ka [sin ^ a » e_,nz (22c) 


where a 0 = In (y) and a n = - for n / 0. 


Then, by (19'), the Fourier coefficients of G(x) is computed as 


ffn — <Jn — ~ 


(23) 


where u n is the convolution sum 


N/ 2-1 

u n = ^ Pm a n—m 


(24) 


m=—N/2 

We note that the convolution sums in (24) require TV multiplications for each u n . 
Thus the total operations for the convolution sums are of order 0{N 2 ). This cost can 
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be reduced considerably to 0(N log 2 N) by the use of a pseudospectral transformation 
method with de-aliasing techniques 7 ’ 8 . For completeness, evaluation of (24) with a 
“padding” de-aliasing technique is given in Appendix B. 


III. SPECTRAL METHOD FOR NEUMANN BOUNDARY CONDITIONS 

We now discuss the Fourier spectral method for the Boundary Integral Equation 

(8b) of the Neumann boundary conditions. Upon substituting the truncated Fourier 

series of the layer distribution function f($) into (8b), we get 
N/2 — 1 N/ 2—1 j. / q2 q 


f E /• 

n=-N/ 2 


e inS + 


n=-N/2 


if 1 r, r 2 'f a 2 G . aa\ „ M ] . ,„»■ 

aM r^ m K 


(25) 


where b n are the Fourier coefficients of the specified Neumann boundary condition. 

Again the integral appearing in equation (25) is directly related to the Fourier 
coefficients of and It is easy to find that the Fourier coefficients of §£ are 

the same as those of already given in the previous section as h n . The apparent 
difficulty here is with the second normal derivative of the Green’s function ^ n i^ n ■ It can 
be shown that this function is strongly singular at x = 0 and, indeed, is not integrable 
in the ordinary sense. Fortunately it can also be shown that the integral with the 
second normal derivative can be transformed into one involving tangential derivatives 
with reduced singularity. In particular, we have 10 


JO 


d 2 G 


in& 


dn'dn 
1 d 


add 


2n fl de in6 1 dG 


I fl 

Jo L a 


+ /c 2 !^ • n G e 


ind 


a dd 


(26) 


flfl a M‘ 

represent tangential derivatives on the boundary. 

The right hand side of (26) is now integrable in the sense of Cauchy Principal Value. 
To show this, we only need to note that by the expression of the Green’s function given 


where and - 

add a dd' 
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in (6) we get 


dG_ _ id 


= ~Tai H ° ) ( 2Ka 


. X 

sin — 

2 



. x 

\ d 

. x 

ZKd 

sin — 

) T7 — 

sin — 

\ 

2 

/ dx 

2 


. x 

\\ sin x 

| Lna 

sm — 


- — 

\ 

2 

\J jsin 

fi 


= 7«, ( " (2 

Recalling (20), the asymptotic expression of for small x is found as 

in a: ko, ( | . x|\ / . x|\ sin a: _ , „ 

- a r° l s,n 2 1) Ji v™ sm si) fun + °w 


( 27 ) 


dG 

d6' 


sin 
87 t I sin 


(28) 


where 0(x ) denotes smooth terms of order x and higher. 

The singular first term shown above is integrable in the sense of the Cauchy Prin- 
cipal Value. In fact, we have 


_L [ 2 * 

2tt J 0 


sin x 


sin f 


e mx dx = 


0 when n = 0 

2 i sign(rc) when n / 0 


(29) 


Upon substituting x — 0 — $' and equating the coefficients of e* n6 , equation (25) 
is reduced to 

+ (^^( X ) + * 2 cos (z) £(*) e' nx adx = b n (30) 

in which we have used the fact that, for a circular cylinder, 


n' • n = cos(0 - 9') 


The integral in (30) will now be evaluated through the Fourier coefficients of each term. 
For the first term, the Fourier coefficients of fp- are obtained from the relation 


dG 

d$' 


dG 

dx 


N/ 2-1 

E in 9ne~ tnx 


n=-N/2 


where g n are the Fourier coefficients of G(x) by (13). 


( 31 ) 
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The Fourier series approximation of the second term in the integral of (30) can also 


be found using g n since we have 


N/ 2-1 Nj 2-1 

cos(x) G(x) = cos(x) <7ne _,ni « ,nx 

n——N/2 n=-N/2 


where 


9n 


2 9-»+l 


n = -N/2 


= l Ug n -l + 9n+i) —N/2 + 1 < n. < N/2 — 2 


2^5r-2 


n = N/2 — 1 

Hence equation (30) is reduced to the following algebraic equations 
y/» + 2 naf n (-^j9„ + it 2 g„ - u/A„) = 4„ 


(32) 


(33) 


(34) 


for —N/2 < n < JV/2 - 1. 

Therefore, the Fourier coefficients of the layer distribution function f(9) for the 
Neumann boundary conditions are obtained explicitly as 

/„ = t^- y (35) 

y + 27 ra \--^9n + K 2 g n - ir]h n ) 

where g n , g n and h n are computed by (23), (33), and (21), respectively. 

We point out, however, that g n as given by (33) and, indeed, h n of (21), are not 
exact for n = — N/2 and N/2 — 1, owing to a truncated series of G(x ) in the compu- 
tation. Whereas it is possible to compute these two coefficients exactly, the resulting 
error in the last two coefficients of /„ is negligible because 6„, in the numerator, de- 
cays exponentially as for smooth boundary conditions. That is, f n for n = —N/2 and 
N/2 — 1 are necessarily negligibly small if N is sufficiently large. For simplicity and 
practicality, (21) and (33) are retained in the numerical calculations. 
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IV. NUMERICAL EXAMPLES 

In this section, numerical results of a plane wave scattering by a circular cylinder 
are presented. The incident wave is assumed to be 


4>i = e 


l fiX 


The scattered wave, <j > , satisfies the Helmholtz equation (1). The boundary conditions 
considered here are the Dirichlet type <j> = and the Neumann type 

The solutions for the scattered field are obtained by the layer representation (2) as 

(^-<vG^mad« 

n / 2-1 2 

= a fnj 


r2 V^-<, G ) 


\ dn 


de 


The above integral can be easily evaluated directly using FFT, since the Green’s func- 
tion has no singularity for points lying outside of the boundary. The details are omitted 
here. 

For plane incident waves, exact solution is given by infinite series of the Bessel and 
Hankel functions®. Our purpose here is to demonstrate the exponential rate of conver- 
gence of the numerical solutions. We emphasize again that the numerical formulation 
applies to any form of the incident waves. Due to its simplicity, a sample FORTRAN 
program is listed in Appendix C. 

In numerical calculations, the radius of the cylinder, a, is taken to be 1 and also 
V = 1* Computations for Ka = 1, 10 and 100 have been carried out. In Tables I to 
IV, numerical values of the layer distribution function f(0) and the scattered function 
(f> at far field are given for selected points in space. Exact values at far field are also 
shown in the tables. Clearly as the number of Fourier collocation points increases, the 
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numerical solution converges exponentially fast. Significant improvements in accuracy 
are observed with relatively small increase of the number of datum points. This is 
often true for spectral methods in general. The error decreases dramatically when the 
number of points is large enough to resolve the basic features of the solution. 

The corresponding layer distribution function f{9 ) is plotted in Figures 1 to 6 
for the Dirichlet and Neumann boundary conditions for Ka — 1, 10 and 100. These 
graphs demonstrate again the remarkable accuracy of the Fourier spectral methods 
with relatively small number of datum points. 

Far field scattered intensities, computed as \r\<f> 2 , are plotted in Figure 8 and 9 for 
the Dirichlet and Neumann boundary conditions, respectively. 


V. CONCLUSIONS 

A fast numerical solution of wave scattering by a circular cylinder has been pre- 
sented. It is shown that by properly removing the non-smoothness of the integral ker- 
nels of the Boundary Integral Equations, Spectrally accurate numerical solutions are 
obtained. The numerical error decays exponentially as the number of datum points 
increase. This implies that the present method requires significantly fewer points for 
achieving a given accuracy in comparison with previous numerical approaches. The 
present method is also easy to implement. 

Moreover, the combined single and double-layer formulation of the Boundary Inte- 
gral Equations ensures the uniqueness of the numerical solution for all frequencies. It 
is shown that the hypersingularity of the Boundary Integral Equations can be handled 
easily in the spectral method. 
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APPENDIX 


A. Exact expressions of g n and h n 

In this appendix we derive the exact analytic expressions for the Fourier coefficients 
of G{x) and §£(x). 

It can be shown that, e.g. by (7.2.51) of ref. 6, 

OO 

tf 0 (1) (2/ca|sin||) = £ J m (Ka)e~ imx (ill) 

m=—oo 

Hence 

1 f 2 * 

9n = — G{x)e' nx dx 
27T Jq 

= kSo ( 2na | sin ||) e '“ ix 

Moreover, for n ^ 0, using integration by part and (Al), 



where use has been made of the formula 9 

zH \\ z )\ = zH o\ z ) 


17 


For n = 0, further calculations show that 


K 2 ai 


A 0 = ^|H< 1) MJ 2 (-co)-/fJ ,) 

K“a . . 1 


16 

,2, 


(«a)J 0 (/ca) - ^H^\Ka)J\{Ka) 


= —\92 - go ) - -91 


B. Evaluation of convolution sums 

An algorithm of computing convolution sums u n with 0{N log 2 N ) operations is 

shown below 8 . 

Let M > 3 N and 

= 2nj/M , j = 0, 1,2, ...M - 1 

Compute the following using FFT for j = 0, 1,2, ...M — 1 . 

M/2-l 


where 


A,= E 


flnje 


-imtj 


m=— M/2 
M/2-1 

p,= E P"-'"'"* 

m=— M/2 


a ra 


— 


—AT < m < JV - 1 


other 


Pm — 


Pm —N/2 < m < N/2 - 1 


other 


and form the product 

Uj — AjPj 

Then the convolution sum u„ is the (backward) FFT of Uj as follows 


M— 1 




j = 0 
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for —N/2 <n<N/ 2-1. 


C. FORTRAN program 

A FORTRAN program of implementing the Fourier spectral method is listed below, 
routines cftti, cfttf and cfttb denote initialization, forward and backward FFT 
transforms respectively.) 

program circle 

c n : number of points; isoft=i : Dirichlet B.C.; isoft=0 : Neumann B.C. 

C*********^********************************************^^^^^^^^^ 

parameter (n= 32 ,ak=i 0 . 0 ,isoit=i,eta= 1 . 0 ,nl=n-l,nhalf=n/ 2 ,m= 3 *n, 

> rn=i loat (n) , pi=3 . 14159265358979324 , euler=0 . 57721666490153286) 
complex b(0:nl) ,ln(0:nl) ,gbar(0:nl) ,gn(0:nl) ,hn(0:nl) ,p(0:nl) , 

> gtilde(0:nl) ,am(0:m-l) ,pm(0:m-l) ,wsave (2000), wsave2(2000),ei, phi 
ei=(0. 0,1.0) 

call cffti(n,wsave) 
call getbc(n,ak,b,ei,pi) 
call cfftf (n,b,wsave) 
do 10 j =0 , n-1 

tmp= 2 . 0*ak*abs (sin(pi*float ( j )/rn) ) 
if(j.eq.O) then 
gbar (0)=-euler/2. 0/pi+ei/4 , 0 

p(0) =1 . 0 

else 

gbar( j )=ei/4. 0*(besj0(tmp)+ei*besy0(tmp) ) 

> +0.5*alog(tmp/2.0)*besj0(tmp)/pi 
p(j )=bes jO(tmp) 

endif 

10 continue 

call cf f tb (n , gbar , wsave) 
call cfftb(n,p, wsave) 
am(0)=alog(ak/2.0) 
am(2*n)=-l . 0/2 . 0/rn 
do 21 i=l ,n-l 

am(i)=-i. 0/2. 0/f loat (i) 

21 am(2*n+i) =1 . 0/2 . 0/f loat (i-n) 
do 22 i=0,nhalf-l 

pra(i)=p(i) 

22 pm(5*nhalf +i) =p(nhalf +i) 
call cffti(m,wsave2) 

call cfftf (m,am,wsave2) 
call cfftf (m, pm, wsave2) 
do 23 j=0,m-l 

23 pm( j )=am( j ) *pm( j ) 
call cfftb(m,pra,wsave2) 
do 31 i=0,nhalf-l 
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gn(i)=gbar(i)-0 . 5*pn(i) /float (m)/pi 

31 gn(nhalf +i) =gbar (nhalf +i)-0. 6*pm(5*nhalf+i)/f loat (m)/pi 
hn(0)=ak**2/4.0*(gn(2)-gn(0))-gn(l) 

hn(nhalf-l)=ak**2/4. 0/float (nhalf-l)*gn(nhalf-2) 

hn (nhalf ) =ak**2/4 .0/float (nhalf ) *gn (nhalf + 1 ) 
hn(n-l)=ak**2/4.0*(gn(0)-gn(n-2)) 
gtilde(0)=0. 5*(gn(l)+gn(n-i) ) 
gtilde (nhalf -1 ) =0 . 5*gn(nhalf-2) 
gtilde (nhalf ) =0 . 5*gn(nhalf +1) 
gtilde(n-i)=0 . 5*(gn(0)+gn(n-2) ) 
do 32 i=l ,n-2 
itrue=i 

if (i.ge. nhalf ) itrue=i-n 
if (i.eq.nhalf-i. or. i.eq. nhalf) go to 32 
hn(i)=-aJc**2/4. 0/float (itrue)*(gn(i+l)-gn(i-i) ) 
gtilde ( i) =0 . 5* (gn(i-l)+gn(i+l) ) 

32 continue 

do 40 i=0 ,n-l 

if (isof t . eq. 1) then 

fn(i)=b(i)/(0.S*rn+2.0*pi*(hn(i)-ei*eta*gn(i))) 

else 

itrue=i 

if (i.ge. nhalf ) itrue=i-n 

fn(i)=b(i)/(0 . S*ei*eta*rn+2.0*pi* (-float (itrue)**2*gn(i) 
> +ak**2*gtilde(i)-ei*eta*hn(i) ) ) 

endif 

40 continue 

c *********************************************************** 

c The following is to find phi at far field r=r0 

C *********************************************************** 
r0=10 . 0 
npoint=4 

do 70 ii=i ,npoint 

s j=2 .0*pi*f loat (ii-l)/f loat (npoint) 
do 71 j=0,n-l 

theta=2 .0*pi*float( j)/rn 
r j =sqrt ( 1 . 0+r0*r0-2 . 0*r0*cos (theta-s j ) ) 
dj=i . 0-r0*cos(theta-sj) 
trap=ak*r j 

gn( j )=ei/4.0*(besj0(tmp)+ei*besy0(tmp)) 

71 hn(j)=-ei*ak/4.0*(besjl(tmp)+ei*besyl(tmp))*dj/rj 

call cf ftb(n, gn, wsave) 
call cfftb(n,hn, wsave) 
phi=0 . 0 
do 72 i=0,n-l 


72 

phi=phi+2. 0*pi*fn(i)*(hn(i)-ei*eta*gn(i) )/m 

70 

write(3 , 100) rO , sj ,phi , cabs (phi) 

phi= * , 3el7 . 10) 

100 

f orxaat ( ’ r0 => ,el8 .6 , ’ theta-* ,615.6/* 

999 

stop 


c 

end 



subroutine getbc(n,ak,b,ei,pi) 
complex ei,b(0:n-l) ,tmp 
do 10 j=0 ,n-l 
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tmp=ei*ak*cos(2 . 0*pi*f loat ( j )/fioat (n) ) 
10 b( j )=-cexp(tmp) 

return 
end 
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TABLE I 

Values of the layer distribution function f(9) at selected points 
on the boundary. Dirichlet boundary condition. 


N 

9 = 0 ° 

9 = 90° 

9 = 180° 

Error 


na = 1 

4 

1.101447573 

1.102982967 

1.124378820 

10" 2 

8 

1.113205176 

1.095419894 

1.146430615 

10" 3 

16 

1.112753432 

1.094877536 

1.145739275 

10“ 8 

24 

1.112753420 

1.094877525 

1.145739263 

10- 12 


na = 10 

24 

4.590213453 

6.904710445 

5.180354736 

10“ 2 

32 

4.546357630 

6.901732036 

5.132718905 

10 -3 

48 

4.545461066 

6.901500667 

5.132515158 

10" 8 

56 

4.545461055 

6.901500659 

5.132515156 

10- 12 


na = 100 

224 

20.64255659 

6.841653547 

18.93255934 

10 -3 

256 

20.64325731 

6.842244857 

18.93221646 

10“ 9 

512 

20.64325733 

6.842244863 

18.93221644 

10- 12 
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TABLE II 

Values of the scattered function at selected points 
at far field r = 10a. Dirichlet boundary condition. 


N 

II 

O 

o 

II 

to 

o 

o 

II 

00 

o 

Q 

Error 


Ka = 1 

4 

0.4146449903 

0.2787718545 

0.1852248716 

10" 2 

8 

0.4224209076 

0.2612785029 

0.2551151985 

10“ 4 

16 

0.4224153154 

0.2613031445 

0.2552183381 

10-1° 

Exact 

0.4224153154 

0.2613031445 

0.2552183381 



na = 10 


24 

0.8255952003 

0.1969679200 

0.1864749710 

10~ 2 

32 

0.8285176644 

0.1953580665 

0.2300067055 

10" 4 

48 

0.8285110664 

0.1953543814 

0.2300939707 

lO-io 

Exact 

0.8285110664 

0.1953543814 

0.2300939707 


Ka = 100 


224 

0.8562228283 

0.1881301853 

0.2295232548 

10 -3 

256 

0.8562289911 

0.1881326409 

0.2294229274 

10-1° 

Exact 

0.8562289911 

0.1881326409 

0.2294229274 
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TABLE III 

Values of the layer distribution function /($) &t selected points 
on the boundary. Neumann boundary condition. 


N 

0 = 0 ° 

0 = 90° 

o 

O 

00 

II 

Error 


Kd = 1 


4 

1.035182633 

0.3028073027 

0.8616587030 

10" 1 

8 

1.200134116 

0.3972281648 

0.8518411247 

IO" 2 

16 

1.199187560 

0.3963806796 

0.8495643896 

io- 7 

24 

1.199187560 

0.3963806589 

0.8495643587 

IO" 12 


a 

II 

I—* 

o 

24 

0.6004486353 

0.4814454225 

1.362228889 

io- 1 

32 

0.6274625969 

0.6575899642 

1.577833267 

10" 2 

48 

0.6302381163 

0.6567081358 

1.460119442 

10- 7 

56 

0.6302381517 

0.6567081198 

1.460119455 

IO" 12 

k a = 100 


224 

0.2185547081 

1.282490390 

2.054272775 

IO -2 

256 

0.2157948725 

1.283008634 

2.057912965 

10- 7 

512 

0.2157947803 

1.283008643 

2.057913072 

10- 12 


24 

























































* TABLE IV 

Values of the scattered function <f> at selected points 
at far field r = 10a. Neumann boundary condition. 


N 

II 

o 

o 

II 

<0 

o 

o 

— 

II 

00 

o 

© 

Error 


Ka = 1 

4 

0.1583300606 

0.1690204144 

0.1619964200 

lO -1 

8 

0.1732916160 

0.1563414831 

0.2312523394 

10“ 4 

16 

0.1733358919 

0.1563260243 

0.2313583724 

io - 10 

Exact 

0.1733358919 

0.1563260243 

0.2313583724 


o 

11 

24 

0.7679584467 

0.2167643382 

0.1574136977 

io - 1 

32 

0.7740714632 

0.1956069424 

0.2282238894 

10~ 3 

48 

0.7740874173 

0.1955960691 

0.2283394143 

10- 10 

Exact 

0.7740874173 

0.1955960691 

0.2283394143 


Kd - 100 


224 

0.7688015277 

0.1871656432 

0.2295250315 

10 -3 

256 

0.7688018590 

0.1871717295 

0.2293995512 

10-io 

Exact 

0.7688018590 

0.1871717295 

0.2293995512 
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Figure 1 . Layer distribution function f(6) for aco = 1 . Dirichlet boundary condition. 
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Figure 2. Layer distribution function f(8) for na = 10. Dirichlet boundary condition. 
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Figure 3. Layer distribution function f(9) for k a = 100. Dirichlet boundary condition. 
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Figure 4. Layer distribution function f{6) for /ca = 1. Neumann boundary condition. 
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Figure 5. Layer distribution function f(9) for k a = 10. Neumann boundary condition. 
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Figure 7. 




Directivities of the far filed scattered function, Dirichlet boundary condition. 
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Figure 8. Directivities of the far filed scattered function, Neumann boundary condition. 
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